function Hmus = mus(Nx,Ny,mu1,mu2)
%MUS Generate the on-site term of Hamiltonian
N1=Nx*Ny; % on-lattice
N2=(Nx-1)*(Ny-1); % off-lattice
N=N1+N2;
Hmus=zeros(N);

for jj=1:N1
    Hmus(jj,jj)=mu1;
end
for jj=N1+1:N
    Hmus(jj,jj)=mu2;
end


end

